The applications of solutions to the H-J equation are numerous. The equation arises in the fields of computer vision, image processing, geoscience, and medical imaging and analysis. For example in computer vision, the shape-from-shading problem, which infers 3D surface shape from the intensity values in 2D image, can be modeled and solved with the Eikonal equation, which is a special form of the H-J equation. Extracting the medial axis or skeleton of the shape can be done by analyzing solutions of the H-J equation with the boundaries specified at the shape contour.
Solutions to the H-J equation have been proposed for noise removal, feature detection and segmentation. In physics, the H-J equation arises in models of wavefront propagation. For instance, the calculation of the travel times of the optimal trajectories of seismic waves is a critical process for seismic tomography. Several methods based on the H-J equation have also recently been introduced as a means for describing connectivity in white matter in medical image analysis.
The Hamilton-Jacobi partial differential equations (PDEs), are given by the following equation (1):H(∇u,x)=1, ∀xεΩ  (1)where Ω is a domain in Rn, u(x) is solution, which can be considered as a travel time or distance from the boundary conditions. Of particular interest is the special formH(∇u,x)=√{square root over ((∇u)M(∇u)T)}{square root over ((∇u)M(∇u)T)}  (2)where M is the speed tensor matrix defined on Ω. We use the Hamiltonian equation defined below for our model equation (3):
                                          H            ⁡                          (                              p                ,                q                ,                r                            )                                =                                                    ap                2                            +                              dq                2                            +                              fr                2                            +                              2                ⁢                                  (                                      bpq                    +                    cpr                    +                    eqr                                    )                                                                    ⁢                                  ⁢                              M            =                          [                                                                    a                                                        b                                                        c                                                                                        b                                                        d                                                        e                                                                                        c                                                        e                                                        f                                                              ]                                ,                      p            =                                          ∂                H                                            ∂                x                                              ,                      q            =                                          ∂                H                                            ∂                y                                              ,                      r            =                                          ∂                H                                            ∂                z                                                                        (        3        )            where p, q, and r are partial derivatives of ui at x along x, y, and z axis, and a, b, c, d, e, and f are upper triangular elements of the matrix M. Equation 1 becomes the Eikonal equation when M is an identity matrix times a scalar function, f( ) which is often called the speed function.
A number of different numerical strategies have been proposed to efficiently solve the H-J equation. These methods can generally be classified into two groups. One group is a class of iterative methods based on a fixed-point update using Jacobi or Gauss-Seidel schemes. Other early work solves the Eikonal equation, a special case of H-J equation, by updating the solutions of the grid using a pre-defined updating order and Godunov upwind Hamiltonian until they converge. This method is simple to implement and produces viscosity solutions, but involves many iterations to converge, and for the worst case situation, complexity can approach the order of O(N2) where N is the number of data elements to be processed. A Fast Sweeping method has also been proposed, which uses a Gauss-Seidel updating order for fast convergence. The Fast Sweeping method has a computational complexity on the order of O(kN) where N is the number of elements to be processed and k depends on the complexity of the speed function. The Fast Sweeping method and a Godunov upwind discretization of the class of convex Hamiltonians can be employed to solve anisotropic H-J equations. Another interpretation of Hamiltonians has been introduced based on the Legendre transformation, which appears to be a Godunov Hamiltonian. This method employs the Lax-Friedrichs Hamiltonian for arbitrary static H-J equations. The proposed method is simple to implement and can be used widely on both convex and non-convex H-J equations, but it involves many more iterations than the Godunov Hamiltonian and the solution shows excessive diffusion due to the nature of the scheme. In general, the iterative methods are slow to converge and are not suitable for interactive applications.
The second group of H-J solvers is based on adaptive updating schemes and sorting data structures. A Dijkstra-type shortest path method has been used to solve convex H-J equations, which are generally referred to as the Fast Marching methods. The main idea behind this method is that solutions for a convex Hamiltonian depend only on the upwind neighbors along the characteristics, so the causality relationship can be determined uniquely and the correct solutions can be computed by only a single pass update. The complexity of the Fast Marching method is O(NlogN), which is the best possible (optimal) asymptotic limit for the worst-possible input data (i.e. worst-case optimal). In this algorithm the running time is slightly affected by the complexity of the speed function. However, for a class of general H-J equations, tracing the characteristics can cause expensive searching among a wider range of neighborhoods of nodes than solving equations using an iterative numerical method. In addition, the method uses a global sorting data structure, e.g., a heap, and therefore the parallelization is not straightforward.